In [1]:
# This changes the current directory to the base saga directory - make sure to run this first!
# This is necessary to be able to import the py files and use the right directories,
# while keeping all the notebooks in their own directory.
import os
import sys
from cStringIO import StringIO
import numpy as np
import requests
if 'saga_base_dir' not in locals():
saga_base_dir = os.path.abspath('..')
if saga_base_dir not in sys.path:
os.chdir(saga_base_dir)
In [2]:
for module in ['hosts', 'targeting', 'mmthecto']:
if module in globals():
reload(globals()[module])
else:
globals()[module] = __import__(module)
g = targeting.get_gama() #re-caches the gama catalo
In [3]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
In [4]:
from astropy import units as u
from astropy.table import Table
from astropy import table
from astropy.coordinates import SkyCoord
In [5]:
# this version includes those both with *and* without incomplete data
mltab = Table.read('catalogs/SAGA.objid.noclean.PROBS_WISE_PROBS.sept28.fits.gz')
mltab
Out[5]:
In [6]:
allspec = Table.read(os.environ['HOME']+'/Dropbox/SAGA/data/allspectaken_v4.fits')
allsats = allspec[(allspec['SATS']==1)&(allspec['REMOVE']==-1)]
allsats
Out[6]:
In [7]:
# info on all the hosts to decide which to investigate further
nsaid2name = {}
for row in allsats:
nsaid2name[row['HOST_NSAID']] = row['HOST_SAGA_NAME']
Table(rows=[(i,
nsaid2name[i],
np.sum(allsats['HOST_NSAID']==i),
np.sum((allsats['HOST_NSAID']==i)&(allsats['MASKNAME']!='sdss')))
for i in np.unique(allsats['HOST_NSAID'])],
names=['NSAID', 'host name', 'n_sats', 'n_sats_not_sdss']).pprint(max_lines=1000)
In [26]:
ody = hosts.odyssey
In [12]:
odycat = Table.read('catalogs/base_sql_nsa147100.fits.gz')
odycat
Out[12]:
In [15]:
odysats = allsats[allsats['HOST_NSAID']==147100]
odysats
Out[15]:
In [20]:
np.all(np.in1d(odysats['objID'], mltab['OBJID_1']))
Out[20]:
OK, all sats are represented, at least
In [22]:
np.sum(np.in1d(odycat['objID'], mltab['OBJID_1'])), np.sum(np.in1d(odycat['objID'], mltab['OBJID_2'])), len(odycat)
Out[22]:
In [25]:
np.sum(np.in1d(mltab['OBJID_1'], odycat['objID'])), np.sum(np.in1d(mltab['OBJID_2'], odycat['objID']))
Out[25]:
In [36]:
mlclosemsk = SkyCoord(mltab['RA'], mltab['DEC'], unit=u.deg).separation(ody.coords) < 2*u.deg
np.sum(np.in1d(mltab[mlclosemsk]['OBJID_1'], odycat['objID'])), np.sum(np.in1d(mltab[mlclosemsk]['OBJID_2'], odycat['objID'])), np.sum(mlclosemsk)
Out[36]:
Lets just forge ahead assuming OBJID_1 is the thing to match against
In [52]:
mltab['objID'] = mltab['OBJID_1']
joinedcat = table.join(odycat, mltab)
joinedsats = table.join(odysats, mltab)
del mltab['objID']
joinedcat
Out[52]:
In [344]:
def find_similar_objects(baseobj, catalog, propertytolerances):
"""
minus sign on propertytolerances means *multiplicative* tolerance (>1). Can also
be a function that takes the base value and comparison vals and gives a mask
"""
msk = np.ones(len(catalog), dtype=bool)
for prop, tol in propertytolerances.items():
if '-' in prop:
m1, m2 = prop.split('-')
compval = catalog[m1] - catalog[m2]
baseval = baseobj[m1] - baseobj[m2]
else:
compval = catalog[prop]
baseval = baseobj[prop]
if callable(tol):
msk = msk & tol(baseval, compval)
elif tol >=0:
msk = msk & (np.abs(compval - baseval) < tol)
else:
msk = msk & ((baseval/-tol)<compval) & (compval<(baseval*-tol))
return catalog[msk]
def show_similar_objects(baseobj, catalog, propertytolerances, nameinfo=None, sorton=None):
simcat = find_similar_objects(baseobj, catalog, propertytolerances)
ras = np.insert(simcat['ra'], 0, baseobj['ra'])
decs = np.insert(simcat['dec'], 0, baseobj['dec'])
if nameinfo is None:
nameinfos = []
else:
nameinfos = [nameinfo] if isinstance(nameinfo, basestring) else nameinfo
names = []
if sorton is not None:
if sorton.startswith('-'):
sorton = sorton[1:]
reversesort = True
else:
reversesort = False
sortvals = []
enum = list(enumerate(simcat))
enum.insert(0, (-1, baseobj))
for i, row in enum:
pieces = ['Realobj'] if i==-1 else ['Fakeobj' + str(i)]
for nameinfo in nameinfos:
label = 'p' if nameinfo.startswith('PROBABILITY') else nameinfo
pieces.append('{}={:.3f}'.format(label, row[nameinfo]))
names.append('_'.join(pieces))
if sorton is not None:
sortvals.append(row[sorton])
if sortvals:
sorti = np.argsort(sortvals)
if reversesort:
sorti = sorti[::-1]
ras = ras[sorti]
decs = decs[sorti]
names = np.array(names)[sorti]
targeting.sampled_imagelist(ras, decs, names=names, n=np.inf, posttoimglist=0.5)
return simcat
In [346]:
tols = {'r':1, 'g-r':.2, 'petroR50_r':-3.5,
'spec_z':lambda objval, catval:(np.abs(catval-objval)/objval > .5) & (catval>-1),
}#'PROBABILITY_CLASS_1': lambda objval, catval: catval < (objval*.5)}
nonsatcats = []
for sat in joinedsats:
nonsatcats.append(show_similar_objects(sat, joinedcat, tols, ['spec_z', 'PROBABILITY_CLASS_1'],
sorton='-PROBABILITY_CLASS_1'))
In [347]:
realsats_idxs = [2, 5, 7]
nonsat_idxs = [[1, 14], [12, 11], [1, 33]]
In [350]:
#convert those to rows
realsats_rows = joinedsats[np.array(realsats_idxs)]
nonsat_rows = [[nonsatcats[ri][li] for li in lis] for ri, lis in zip(realsats_idxs, nonsat_idxs)]
In [351]:
def row_to_cutout_img(row, size=(400, 400), scale=0.1):
baseurl = 'http://skyservice.pha.jhu.edu/DR12/ImgCutout/getjpeg.aspx'
res = requests.get(baseurl, dict(ra=row['ra'], dec=row['dec'], scale=scale, width=size[0], height=size[1]))
sio = StringIO(res.content)
sio.seek(0)
return mpimg.imread(sio, 'jpeg')
In [359]:
fig, spls = plt.subplots(len(nonsat_idxs[0])+1, len(realsats_idxs), figsize=(10, 10))
for ax, row in zip(spls[0], realsats_rows):
imdata = row_to_cutout_img(row)
ax.imshow(imdata)
ax.set_xticks([])
ax.set_yticks([])
for imcol in range(len(nonsat_idxs[0])+1):
for ax, row in zip(spls[1:, imcol], nonsat_rows[imcol]):
imdata = row_to_cutout_img(row)
ax.imshow(imdata)
ax.set_xticks([])
ax.set_yticks([])
spls[0, 0].set_ylabel('Real Satellites', fontsize=20)
for spl in spls[1:, 0]:
spl.set_ylabel('Non-Satellites', fontsize=20)
fig.tight_layout()
plt.savefig('satnonsatcomp2.png')
plt.savefig('satnonsatcomp2.pdf')
In [360]:
print 'Real sats:'
for row in realsats_rows:
print row['ra'], row['dec'], row['spec_z'], row['PROBABILITY_CLASS_1']
In [365]:
print 'Non-sats:'
for i,rows in enumerate(nonsat_rows):
print 'For realsat #', i+1
for row in rows:
print row['ra'], row['dec'], row['spec_z'], row['PROBABILITY_CLASS_1']
In [295]:
fig, spls = plt.subplots(len(lowprob_idxs[0])+1, len(realsats_idxs), figsize=(7, 12))
for ax, row in zip(spls[0], realsats_rows):
imdata = row_to_cutout_img(row)
ax.imshow(imdata)
ax.set_xticks([])
ax.set_yticks([])
for imcol in [0, 1]:
for ax, row in zip(spls[2:, imcol], lowprob_rows[imcol]):
imdata = row_to_cutout_img(row)
ax.imshow(imdata)
ax.set_xticks([])
ax.set_yticks([])
spls[0, 0].set_ylabel('Real\nSatellites', fontsize=20)
spls[1, 0].set_ylabel('High-Prob\nNon-Satellites', fontsize=20)
fig.tight_layout()
spls[3, 0].set_ylabel('{0}Low-Prob\n{0}Non-Satellites'.format(' '*36), fontsize=20)
plt.savefig('satnonsatcomp.png')
plt.savefig('satnonsatcomp.pdf')
In [11]:
cat150887 = Table.read('catalogs/base_sql_nsa150887.fits.gz')
In [15]:
allspec5 = Table.read(os.environ['HOME']+'/Dropbox/SAGA/data/allspectaken_v5.fits.gz')
allspecin150887 = allspec5[allspec5['HOST_NSAID']==150887]
In [17]:
mltab['objID'] = mltab['OBJID_1']
joinedcat = table.join(cat150887, mltab)
joinedspec = table.join(allspecin150887, mltab)
del mltab['objID']
In [23]:
goodspec = joinedspec[joinedspec['ZQUALITY']>3]
goodspec
Out[23]:
In [34]:
dict([(utn, np.sum(utn==goodspec['TELNAME'])) for utn in np.unique(goodspec['TELNAME'])])
Out[34]:
In [41]:
bins = np.linspace(0,1,101)
for utn in np.unique(goodspec['TELNAME']):
plt.hist(goodspec['PROBABILITY_CLASS_1'][utn==goodspec['TELNAME']], histtype='step', label=utn, bins=bins, log=True)
plt.legend(loc=0)
Out[41]:
In [43]:
np.sum(joinedcat['PROBABILITY_CLASS_1']>0.01), np.sum(joinedspec['PROBABILITY_CLASS_1']>0.01)
Out[43]: